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The standard star formation scenario is based on spherically symmetric ra- 
diative hydrodynamics calculations. In this paper, we focus on the evolution of 
the rotating non-spherical cloud through the first core and second collapse phases 
[ 10~ 13 < p < 10 3 g cm -3 using axisymmetric numerical calculations. The rotat- 

ing first core evolves under the influence of mass and angular momentum accretion 
from the infalling envelope. The structure of the rotating first core is character- 

Q\ [ ized by the angular momentum J core and mass M core , both of which increase with 

time j CO re, M corc > 0. The evolution path of a rotating first core is considered as 
a sequence of equilibrium solutions with a constant J COTe /M^ ore . Such evolution 

^ | paths are in good agreement with the results of three-dimensional numerical sim- 

ulations. Only the slowly rotating first core with J coro /M c 2 orc < 0.015G/(v / 2ci SO ) 
begins the second collapse directly. On the other hand, the cloud with a rotation 
rate of J cove /M^ om > 0.015G/(v2ci SO ) does not increase in the central density 
beyond a certain maximum value. Although the first core grows in mass by ac- 
cretion, this first core does not begin the second collapse directly. Such a first 
core has a very long time-scale and seems to suffer from a non-axisymmetric 
evolution, such as formation of massive spiral arms, deformation into a bar, or 
fragmentation. In the usual case, the time-scale of a rotating first core becomes 
several thousand years or more. So we expect that at least several percent of the 
prestellar cores contain a first core. The rotating first core has a small average 
luminosity of L core = 0.003 — 0.03(M core /M re f) L . However, we expect that the 
rotating first core repeats temporal increases in the luminosity during its first 
core phase. 

Subject headings: theory — hydrodynamics — methods: numerical — accretion 
disks — stars: formation 
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1. INTRODUCTION 

It is known that stars are formed through the gravitational collapse of the molecular 
cloud core with a density rise of about fifteen orders of magnitude. A dense pre-stellar core 
has a central density of n ~ 10 5 cm -3 (Motte, Andre, & Neri 1998; Kandri et al. 2005) and 
shows statistically significant inward motions (see, e.g., Zhou et al. 1993; Lee, Myers, & 
Tafalla 2001). It is believed that the prestellar core forms a protostar which is observed as 
a class source in a time-scale of ~ several times of 10 5 yr by statistical studies (Onishi et 
al. 1998). 

However, their large optical depth, lack of tracer, and short time-scale prevent us from 
observing objects between the pre-stellar core (n ~ 10 5 cm~ 3 ) and the protostar (class 
object) phases. The protostar comes to be observed as a class source which had already 
accreted significant mass from the envelope by that time. This missing link is one of the 
most important targets of next-generation observational instruments. Theoretically, it is 
necessary to make clear the process of stellar core formation from prestellar core. This gives 
the initial condition of the main accretion phase afterwards. 

The gstandardh star formation scenario from a prestellar core to a protostar is based 
on the numerical simulations of spherical collapsing clouds core. Larson (1969) and Appen- 
zeller & Tscharnuter (1974) calculated the spherical collapse from gas cloud to protostar 
by using the diffusion (or Eddington) approximations and/or the gray approximation of ra- 
diation transfer. Masunaga, Miyama, & Inutsuka (1998) followed the evolution using the 
radiative hydrodynamical simulation and obtained similar results. The collapsing cloud is 
initially optically thin to the thermal emission from dust grains, and the cloud experiences 
isothermal runaway collapse. During this phase, the cloud collapses in a self-similar manner 
(Penston 1969; Larson 1969). After that, the isothermal condition is broken at the density 
of nR2 — 10 10 cm~ 3 and the gas temperature begins to increase. The collapsing cloud forms 
a hydrostatic equilibrium object with a radius of a few AU and mass of ~ O.O1M , which 
is called the gfirst core.h The first core increases its mass by accretion from the isothermal 
envelope and it raises the central density and temperature. The hydrogen molecules begin 
to dissociate into atoms after the central temperature reaches about 2000 K (n H 2 — 10 16 
cm -3 ). Since this reaction extracts the thermal energy from gas, a gsecond collapseh begins. 
After the dissociation is completed (n m ~ 10 20 cur 3 ), a gsecond coreh(or gstellar coreh) is 
born at the center. 

However, interstellar cloud cores show rotation motions. Goodman et al. (1993) observed 
43 cloud cores by NH 3 (J, K) = (1,1) transition (nu2 — 10 4 cm~ 3 ) and found that the angular 
velocity f2 c of these cores is distributed in the range of (0.3 — 4) x 10 -6 yr -1 . Caselli et al. 
(2002) obtained similar results of Vt c = (0.5 - 6) x 10~ 6 yr" 1 with N 2 H+ data (n m ^ 
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10 5 cm -3 ). These rotations correspond to the ratio of the rotational energy to the self- 
gravitational energy of T/| W| ~ 10~ 4 — 0.07. Rotation of the cloud core must have a crucial 
influence on star formation (for a review see Bodenheimer et al. 2000). 

Even today, computational power is not sufficient for multi-dimensional radiative hy- 
drodynamics, and therefore a number of simplifications are supposed in such investigations. 
Isothermal gas is a good approximation for the gas of nH2 < 10 10 cm~ 3 as well as the spherical 
case. Cloud cores near hydrostatic equilibrium are hard to fragment during the isothermal 
runaway collapse phase (Boss 1987; Boss & Myhill 1995; Miyama 1992; Tsuribe & Inutsuka 
1999). During this phase, the central region of the collapsing cloud spins up and it forms a 
dynamically collapsing disk (Norman, Wilson, & Barton 1980; Miyama, Hayashi, & Narita 
1984). Matsumoto, Hanawa, & Nakamura (1997) have shown that such a collapsing isother- 
mal cloud tends to have a central angular velocity of Q c ~ (0.2 — 0.3) /ts in the runaway 
collapse phase, where ts is the free-fall time scale at the central density. 

After the isothermal condition is broken, many authors follow the subsequent evolution 
by assuming the polytropic equation of state based on the results of spherical cloud collapse 
(Bonnell & Bate 1994a; Bate 1998) or radiative transfer in diffusion approximation (Boss 
1986; Boss & Myhill 1995). Matsumoto & Hanawa (2003, hereafter MH03) followed the 
evolution using the nested grid method with a polytropic equation of state. The first core 
accretes mass and angular momentum from the infalling envelope and thus the radius and 
mass of the first core are much larger than the non-rotating spherical model. The radius 
and mass of the rotating first core are equal to a few tens AU and a few times 0.01 M , 
respectively. They found that the first core also increases in T/|W| by mass accretion. The 
polytropic gas disks with T"/|W| > 0.27 suffer from axisymmetric instability (Durisen et 
al. 1986; Pickett, Durisen, & Davis 1996; Imamura, Durisen, & Pickett 2000) under a wide 
range of conditions. Bonnell & Bate (1994a) confirmed that the first core with T/|W| > 0.3 
becomes unstable against the non-axisymmetric instability to both m — 1 and 2 modes 
during the accretion phase and forms massive arms. Subsequently, the first core increases in 
its density by angular momentum transport via gravitational torque or fragmentation into 
binary or multiples. MH03 did an extensive parameter survey of the collapse of rotating 
Bonner-Ebert clouds with a central density of nH2 — 2.6 x 10 4 cm -3 . They found that first 
cores have a wide variety of morphologies depending on the initial rotation speed, VL c ts and 
that clouds with Q c0 > 0.05/ts fragment into binaries or multiples in the first core phase. 

Thus, the accretion of mass and angular momentum leads to a crucial influence on 
the evolution of the first core and first cores must have a wide variety of evolution paths 
depending on the angular momentum of the infalling matter. However, it is not clear how 
the evolution path of the first core is systematically affected by mass accretion. In this paper, 
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we show the evolution of the first core by using equilibrium solutions comprehensively and 
quantitatively. 

In a non-rotating case, the first core begins to collapse again, and form a stellar core. 
This second collapse is important in understanding star formation, because it presents the 
direct initial condition of the protostar and protostellar disk and it may lead to formation of 
a close binary. However, a lot of obscurities and discussions remain concerning the evolution 
of a rotating cloud. Boss (1989) and Bonnell & Bate (1994b) calculated the second collapse 
phase from a gravitationally unstable rotating first core, assuming a polytropic equation of 
state. They found that the fragmentation occurs through two mechanisms once the collapse 
is halted. They found two typical evolutions: first, a protostellar disk fragments by non- 
axisymmetric m = 1 and 2 modes. Second, a collapsing cloud bounces centrifugally and 
forms a ring at the end of the second collapse phase. This ring results in fragments of several 
pieces. However, they assumed a first core with uniform density and Gaussian density profile 
instead of a rotational equilibrium configuration. This must affect the evolution. 

On the other hand, Bate (1998) obtained a different result. He followed the cloud col- 
lapse from the isothermal phase to the second core formation by the SPH method, assuming 
a polytropic equation of state. In his calculation, the first core forms two massive arms dur- 
ing the first core phase and a large amount of angular momentum is removed from the core. 
This results in a slowly rotating core and forms a stellar core directly, without any fragmen- 
tations, similar to the non-rotating spherical collapse. However, the number of models in 
his letter is so restricted, it is difficult to understand how widely these results are realized. 
We examine hydrodynamical evolution of the second collapse from the end of the first core 
phase, which is obtained by our quasi-static calculation. From extensive parameter study of 
rotation parameters we show the evolution is classified into three types according to rotation 
and density. 

The phases of this paper are as follows: the model of the rotating first core is described 
in §2, and we show the quasi-static evolutional path of the first core with mass accretion 
from the envelope in §3. In §4, we explore the second collapse more comprehensively with 
hydrodynamical simulations. We discusses the results in §5. 

2. ROTATING FIRST CORE 

2.1. Basic Equations and Numerical Method 

We calculate the equilibrium state of the rotating first core with various central densities 
and rotation speeds in order to obtain its evolution path as a series of solutions. The basic 
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equations are as follows, 

(pril 2 , 0, 0) - VP - pVtp = 0, (1) 

and 

A*0 = AnGp (2) 

where v, p, ip, and G denote the velocity, thermal pressure, gravitational potential, and 
gravitational constant, respectively. All the variables are assumed to be functions of r and 
z, and independent of ip in the cylindrical coordinates (r,ip,z). 

For simplicity, we did not solve the energy equation which must include radiative cooling 
and chemical reaction. Instead, we assume a polytropic equation of state, 

P = Kp\ (3) 

where 7 and K are the polytropic index and constant, respectively. To mimic the thermal 
history of the first core, we switched 7 and K at the H 2 dissociation density, pdiss = 5.6 x 1CT 8 
g cm -3 , as follows, 

h K) = { (7/5 ' 5 ' 697 X 1Ql3) 9 < Pdiss ' U) 

1 ' I (1.1, 3.823 x 10 11 ) Pdiss < p , 1 1 

in cgs unit. The temperature is proportional to T oc p 7 " 1 in this polytropic model. The 
transition density and polytropic parameters are taken from the radiative hydrodynamical 
calculations of Masunaga, Miyama, & Inutsuka (1998). 



2.2. Models 

The rotating first core is expressed as an axisymmetric equilibrium of a rotating poly- 
tropic gas cloud with a mass of M core and total angular momentum of J CO re- We assumed 
that the distribution of angular momentum of the rotating first core is the same as one for 
a uniformly rotating, uniform-density sphere, 

where j(M(R)) is the local specific angular momentum at cylindrical radius R from the 
rotational axis, M(R) is the mass within the cylindrical radius of R. 

Density and rotation velocity distribution of the first core are numerically constructed 
by the self-consistent field (SCF) method of Hachisu (1986), which iterates between solutions 
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of the Poisson equation and that of the equilibrium equation among pressure, centrifugal, 
and gravitational forces. We extend the original SCF code to include the effect of polytropic 
index 7 dependent on density (see eq. [4]). The density is assumed to equal to zero outside 
the first core for simplicity. 

The equilibria of rotating first cores are specified by only two model parameters, the 
central density, p c0 , and the total angular momentum, J core , according to Hachisu (1986). 
The mass of the first core M core changes depending on these two model parameters. Figure 1 
illustrates the density distributions of the first cores with the same central density as p c o = 
4pdiss but three different angular momenta. They are typical initial conditions of evolution 
models. Their total angular momenta are J core = 9.99 x 10 49 , 4.18 x 10 49 , and 2.25 x 10 49 
g cm 2 s" 1 and hereafter we call them Models A, B, and C, respectively. The contour curves 
and the gray scale show the logarithm of the density. White dashed curves denote the iso- 
density contour at p = pdiss- Figure 1 shows that all first cores form density plateaus in the 
central part and have concave-shaped envelopes. The first core is larger and more massive 
for models with larger angular momentum. Their masses are M core = 4.97 x 10 31 , 3.45 x 10 31 
and 2.77 x 10 31 g. The thicknesses of their cores are identical 0.4 AU because the cores are 
supported only by the pressure force in the z-direction. 



3. Evolution Path of Rotating First Core with Mass Accretion 

During the first core contraction phase, the mass and the central density of the rotat- 
ing first core increases gradually by means of mass accretion from the infalling isothermal 
envelope. 

Figure 2 illustrates the relation between the central density, p c0 , and the mass of the 
first core, M core . The thick curve denotes a series of non-rotating solutions with J core = 0. 
The non-rotating first core evolves along this curve from left to right as the mass increases 
by accretion from the infalling envelope. When it reaches the maximum mass of ~ 0.01 
Mq at p c o — 2pdi SS , further mass accretion results in a dynamical collapse, since there is no 
hydrostatic equilibrium beyond this mass. The cloud goes into the second collapse phase 
(see Larson 1969). This comes from the dependency of the Jeans mass, Mj oc T 3 / 2 p -1 / 2 oc 
p(37-4)/2_ rj^g mass Q f fl rs t core increases in proportion to p^ 1 in the low density part 
p c0 < 10~ 7 g cm~ 3 and decreases in proportion to p~®' 35 in the high density part p c0 > 10~ 7 
g cm~ 3 . 

In the rotating case, the first core becomes more massive by support of the centrifugal 
force against gravitational force. The evolution path is not only controlled by mass accretion 
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but also angular momentum accretion from the infalling envelope. Thus, there are innumer- 
able evolution paths which have different increases in angular momentum J core . Here, we 
assume the mass accreting first core evolves while maintainin the value of J core /M^ ore . It 
is because many numerical simulations show that the specific angular momentum in the 
central region of the rotating cloud tends to be proportional to the enclosed mass during 
the isothermal collapse (Matsumoto, Hanawa, & Nakamura 1997). We introduce a non- 
dimensional rotation rate, u, as an index of the rotation defined as 



V^2C; so t/ C ore / « \ 

ws — ACT' (6) 



core 



where c iso is the isothermal sound speed. Each solid curve in Figure 2 denotes the series of 
the equilibrium solution with the same ui and denotes the evolution paths of the first core 
with the infalling isothermal envelope. 

Figure 2 indicates that the central density increases with the core mass for p c0 < Pdiss- In 
the case of u < 0.015, the rotating first core has a maximum mass at p c0 = (2— 4)pdiss as well 
as the non-rotating case. These is no hydrostatic equilibrium if the mass is larger than this 
maximum value. Such a first core begins the dynamical collapse, which is called the second 
collapse after the mass exceeds the maximum mass. The solutions with dM/dp c o\ w < are 
unstable, and are indicated with a region colored in red in Figure 2. 

In contrast, the mass of the first core increases monotonically with the central density 
for u) > 0.015 and the unstable region dM/dp c \ w < disappears. Thus, the first core 
cannot begin the second collapse, even if M core > O.O18M at p c o — 4pdiss- This means that 
the maximum mass of the core that experiences the second collapse is 0.018 M Q . The central 
density hardly increases even if the accretion increases the mass. Such a first core dose not 
begin the second collapse unless the core spinning slows down owing to the fragmentation 
or the angular momentum transport. 



4. SECOND COLLAPSE PHASE 

4.1. Numerical Method and Boundary Conditions 

We obtained the quasi-static evolutionary path of the first core caused by the mass 
and angular momentum accretion. The end points of the quasi-static evolution have been 
shown in §3. Now, we calculate the hydrodynamical evolution of the rotating first core by 
numerical simulations. The initial conditions are the rotating first core obtained in §2 and 
§3. We examine the stability of the first core and follow the gravitational collapse during 
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the second collapse phase. We superpose 1% density enhancement on the cloud and assume 
a low-density gas of ~ 10~ 5 p diss outside the first core. 

The basic equations for hydrodynamical calculation are as follows, 

|| + V-(p«) = 0, (7) 

d 

— (pv) + V-(pv®v) + VP + pW> = 0, 8) 
ot 

and 

Atp = AixGp. (9) 

All the variables are assumed to be functions of r and z, and independent of ip in the 
cylindrical coordinates (r, <p, z). For simplicity, we did not solve the energy equation and we 
assume a polytropic relation similar to equation (3) as 

= l(T 13 gcni" 3 , 

5.6 x 1CT 8 gem -3 
10~ 3 g cm 



7 = < 
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where pi SO and p ste i are the transition density from the isothermal to adiabatic equation of 
states and the density after which the hydrogen molecule dissociation is completed. 

The calculations were performed using a two-dimensional axisymmetric Euler code. We 
developed the numerical code to solve the hydrodynamical equations (7)-(9) by modifying 
that of Saigo, Matsumoto, & Hanawa (2000). Our code consists of two parts, the solver of 
the hydrodynamical equations and that of the Poisson equation. The former is based on 
the Flux Difference Splitting (FDS) scheme of Roe (1981). The second order accuracy is 
achieved by the Monotone Upwind Scheme for Conservation Law (MUSCL) method (e.g. 
Hirsch 1990). We solved the Poisson equation with the five-point central difference scheme 
and the multi-grid iteration (see, e.g. Press & Teukolsky 1991). 

We achieve high spatial resolution throughout the computations by rezoning. As the cen- 
tral density increases with the gravitational collapse, the central high-density region shrinks. 
To compute the small-scale structure, we resize the grid cell and restrict the computation 
domain to the central part with the grid number of (n r ,n z ) = (512,512) in the r- and 
^-directions. The size of the initial computational domain is taken to be twice as large as 
the disk radius. 

The first core is assumed to be isolated and we assume no flux across the outer boundary. 
The boundary condition of the gravitational potential is set by the spherical harmonics 
expansion. 
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We introduce a new finer grid so that the resolution satisfies the condition that the 
half-width of the half-maximum density should be resolved with more than 16 grid spacings. 
The size of the finer computational domain is taken to be half that of the coarser one. When 
a new grid is introduced, the density and velocity on the grid points are interpolated by the 
first order bilinear method. The density, velocity and gravitational potential are fixed on 
the outer boundary after rezoning. The Poisson equation is solved on the fine grid with the 
boundary gravitational potential at the rezoning epoch. 

4.2. Typical Evolutions 

Here, we show the dynamical evolution of the rotating clouds with a wide range of 
rotation speeds from the end of the first core to the stellar core. The evolution of the first 
core shown to be is divided into three types according to rotation speed. 

4-2.1. Model A: Direct Stellar Core Formation 

In this subsection, we show the second collapse of the slowly rotating cloud of model 
A. It is a typical example in which the runaway second collapse continues till the stellar 
core formation. The parameters of model A are p c o = 4pdiss = 2.3 x 10~ 7 g cm~ 3 and 
uj = 0.012, which has the maximum allowable mass with the same rotation rate, u (Fig. 
3 a). This model shows evolutions qualitatively similar to those of the non-rotating spherical 
collapse model (Larson 1969; Masunaga, Miyama, & Inutsuka 1998). 

Figure 36 shows the density distribution at the age of t — 4105.3 days after the dy- 
namical contraction begins. The first core begins runaway collapse with the gravitational 
instability and the central density p c is ten times larger than that of the initial condition of 
dynamical evolution, p c0 . In this paper, we defined the collapse in which the central density 
rises by one order of magnitude as the second collapse. In the early stage of second collapse, 
the first core collapses mainly vertically. Note that the dynamical collapse is seen not only 
in the central dissociation region p > p diss but also in the outer part of the first core. 

Further evolution was traced with finer grids. Figure 3c and 3d shows the density 
distributions at the ages of t — 4519.3 days and t = 4531.7 days, respectively. The 
central densities are 100 and 1000 times larger than the initial one, respectively. The gas 
disk evolves self-similarly during the runaway collapse phase. The central high density disk 
becomes thinner and thinner during the runaway collapse. This second collapse phase ends 
and the increase in the central density stops at t — 4546.5 days (Fig. 3e). The white 
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solid curve denotes the stellar core of p c > p stc \- The stellar core bounces outwards on the 
equatorial plane (see Bonnell & Bate 1994b). At the stage of t — 4601.8 days, the density 
of the core has fallen below the stellar one and the core is surrounded by a shock wave with 
a radius of 0.02 AU (Fig. 3f). 

Figure 4a shows the density profiles on the equatorial plane. It is seen each density 
profile consists in a plateau around the center and a power-law envelope and the evolution 
proceeds in a self-similar fashion. During the runaway collapse phase (solid curves), the 
density in the envelope is proportional to r~ 2 ' 68 . This power-law index is steeper than that 
of the spherical collapsing polytropic gas of r _2 ^ 2 ~ 7 ^ = r~ 2 ' 22 (7 = 1.1) (see e.g. Suto & Silk 
1988). It is due to the increase in the ratio of the centrifugal force to the pressure gradient 
force during runaway collapse. 

During calculation, Toomre Q-value inside the central plateau of the collapsing cloud is 
below unity. We also plot the unstable wavelength of Aj in figure 4a which is marked as a 
cross. The unstable length is derived from the Toomre reference analysis for the centrifugally 
supported disk, 

2r 2 

Aj = (11) 

G£(l + v/T^) 

where c s , E, and Q denotes the sound speed, surface density of the disk, and Toomre Q-value, 
respectively. It is equal to or less than the radius of the central high density plateau on each 
profile, so the collapsing disk is stable against the axisymmetric instability. 

Figures 4b and 4c show the infall and rotation velocity profiles on the equatorial plane. 
They have maxima at the edge of the central density plateau and decrease with radius in the 
envelope. This distribution is the same as the self-similar solution of the spherical collapsing 
polytrope (Yahil 1983; Suto & Silk 1988). Figure 4 d shows the infall velocity profiles along 
the z-axis. The maximum infall velocity attains 2.5 x 10 5 cm s" 1 at the edge of the central 
high density region. 



4-2.2. Model B: Formation of Rotationally Supported Core 

Even when the rotation is a little faster than model A, the evolution becomes quite 
different from model A. We show model B as a typical example for runaway collapse models 
with larger angular momenta. The initial condition of model B (Fig. 5 a) is p c o = 4pdiss and 
M corc = 3.45 x 10 31 g (u = 0.014). 

Figures 56 and 5c show the density distribution at t — 4982.9 and 5712.9 days. The 
central density of these stages are 10 and 100 times higher than the initial state. As in model 
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A, the cloud begins its runaway collapse as a whole with the gravitational instability. 

In this model, the second collapse suddenly stops and the cloud forms an equilibrium 
core (or disk) before the second core formation. Figure 5d shows the disk (inside the white 
dotted curve) at the stage of p c = 2.45 x 1CT 5 g cm~ 3 which is much lower than the stellar 
density of p stc i = 1.0 x 10~ 3 g cm~ 3 (see eq. [10]). Here, we define the equilibrium core 
as the area where the infall speed becomes less than 20% of sound speed, which may be 
regarded as the 1.5th core. The core is mainly supported by centrifugal force. The force 
fraction of the centrifugal to the pressure force is 0.65. The radius and mass of the 1.5th 
core is 0.02 AU and 8.8 x 10 30 g, respectively. At the stage of Figure be, the radius and mass 
of the disk increase to 0.18 AU and 1.05 x 10 31 g due to mass accretion, respectively. On 
the other hand, the density of the central region decreases due to the centrifugal bounce. 

Figure 6a shows the density profiles on the equatorial plane, p(r). This shows that the 
evolution proceeds roughly self-similarly as well as model A. During the runaway collapse 
phase (solid curves), the density in the envelope is proportional to r~ 2 - 9 . This profile is 
steeper than that of model A. During the accretion phase (dashed curves), the central density 
decreases to about one-half of the maximum value at t — 5734.8 days due to a bounce of 
the disk. The centrifugally supported disk accretes mass through a shock front which is seen 
at r ~ 0.16. 

Figures Qb and 6c show the radial infall, v r , and rotation velocities, v v , on the equatorial 
plane. The evolution is almost similar to model A. At stage (d), the infall velocity has 
vanished inside the radius of 0.02 AU. It is the radius of the high density plateau. Driven 
by the accretion, the disk radius increases to 0.16 AU at stage (e). 

4.2.3. Model C: Stable First Core 

In faster rotating model C, the rotating first core indicates an oscillation around the 
initial equilibrium state and does not show any indications of the runaway collapse. Model 
C is a typical example of such a stable first core. The parameters are p c o = 4pdiss and 
M core = 4.97 x 10 31 g (see Fig. 2). 

4.3. Classification of Second Collapse Phase 

The initial conditions of hydrodynamical calculations are plotted with crosses, triangles, 
and circles in Figure 2. First cores marked by open crosses collapse dynamically until forma- 
tion of the stellar core as shown in model A. These models are distributed below the curve 
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of to — 0.012 in Figure 2. First cores marked by triangles collapse dynamically but form 
centrifugally supported disks below the stellar density. Rapidly rotating cores indicated by 
crosses do not experience the second collapse phase. This non-collapse model is defined as 
the model in which the central density does not increase to 10 times that of the initial one. 
The boundary between the collapse model and non-collapse model roughly corresponds to 
the boundary which separates the stable dM/dp c > and unstable dM/dp c < regions 
already shown in §3 (red region in Fig. 2). 

5. DISCUSSION 
5.1. T/\W\ 

Figure 7 is the same as Figure 2 but it shows the rotational energy of the first core. The 
color scale shows the ratio of rotational energy to the gravitational energy, The first 

core with larger rotational energy T/|W| can support more mass. A dashed curve denotes 
the critical curve of T/|W| = 0.273 for non-axisymmetric instability. That is, Imamura, 
Durisen, & Pickett (2000) showed that the rotating polytropic gas disk with 7 = 7/5 is 
unstable against the non-axisymmetric instability if > 0.273 (see also, Durisen et al. 

1986; Pickett, Durisen, & Davis 1996) 

The first cores with uo > 0.015 stop increasing central density, and the evolution 
path rises nearly vertically with mass accretion. During the evolution, increases 
monotonically along the path. Finally exceeds the critical value 0.273 and the first 

core enters the non-axisymmetric unstable region. On the other hand, the first core with 
uj < 0.015 evolves to the right in Figure 7 and begins the second collapse directly without 
non-axisymmetric instability. In this case, T"/|W| increases slowly and there is a maximum 
value of 0.2 at most. 



5.2. Comparison with Three-dimensional Calculations 

A number of authors have studied the formation and evolution of the first core from the 
rotating cloud with various initial conditions of rotation and density distribution using three- 
dimensional calculations (Bonnell & Bate 1994a; Bate 1998; Saigo, Hanawa, & Matsumoto 
2002; Matsumoto & Hanawa 2003). They found that the first core shows a wide variety of 
evolutions depending on its initial condition. However, its evolution had not been understood 
comprehensively. Here, we show that the variety of results is understood using the central 
density core mass (p c o — M core ) diagram (Fig. 7). 
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Numerical simulations have show that the first core which is formed from a rotating 
cloud becomes non-axisymmetrically unstable in a short time after the formation of the 
first core. This can be understood if we think that all the clouds in their calculations have 
sufficiently fast rotation of uj > 0.015. 

Bate (1998) calculated the collapse of a uniformly rotating (Q c = 7.6 x 10~ 14 s" 1 ) 
uniform density (p = 10~ 18 g cm" 3 ) sphere using the SPH method with a polytropic equation 
of state. The rotating cloud forms a hydrostatic first core with a mass of M core ~ O.O1M 
at p c ~ 2 x 10" 11 g cm -3 (lower square symbol in Fig. 7). After that, a non-axisymmetric 
unstable m = 2 mode begins to grow at p c ~ 10~ 10 g cm -3 . At that time, the first core 
has M corc ~ 0.04M Q , R COTC ~ 40 AU and T/\W\ ~ 0.34 (upper square symbol in Fig. 7). 
Between these stages, the cloud evolves nearly vertically in Figure 7. The evolution almost 
agrees with our evolution path with uj = 0.06. This value roughly agrees with the initial 
value uj ~ £l c ts = 0.08 at p = 10~ 18 g cm" 3 (See Appendix A). 

Afterward, this first core forms two massive spiral arms and the angular momentum is 
transferred to the outer disk (see Figs. 3 and 5 of Bate 1998). The central region begins the 
second collapse and it continues till the stellar core formation. In our analysis, the first core 
with uj = 0.06 can start the second collapse if J COTC / M core 2 of the core becomes 1/4 or less of 
the value before the formation of spiral arms. Bate (1998) reported that the total angular 
momentum of the central disk within several AU decreases to ~ 20% of the value before the 
formation of spiral arms 1 . 

MH03 calculated the formation and evolution of the first core from the rotating Bonner- 
Ebert sphere for 225 models with different initial rotation speeds, rotation laws, and ampli- 
tudes of density perturbation. They found that a cloud core with non-dimensional angular 
speed of fl c ts < 0.03 at the central density of p c = 1.0 x 10" 19 g cm" 3 does not fragment in 
their calculation. These first cores deform into spirals or a bar due to the non-axisymmetric 
instability just after the first core formation. 

In Figure 7, we plot the evolution path of their High-Resolution-Model with Q c tg = 
0.0276 initial at the central density of p c = 1.0 x 10~ 19 g cm -3 (hereafter we refer it as 
'HRM003') by the dashed curve (MH03). The central density increases gradually due to 
mass accretion for p < 10" 9 g cm" 3 . The evolution of HRM003 is almost along our evolution 
path of uj = 0.04. After the evolution path goes over our critical curves for non-axisymmetric 
instability, massive spiral arms appears in the first core (Fig. 3d of MH03). Soon after, the 
central region of the spiral arms shrinks into a compact centrifugally supported core and the 



1 The exact change of the rotation rate oj angular has not been obtained from Bate (1998) owing to a lack 
of information about the mass within several AU. 
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central density increases by an order of magnitude (Fig. 3e of MH03). Thus, their result 
of HRM003 agrees well with our evolution path of uj = 0.04 and it also agrees with our 
stability analysis for non-axisymmetric instability. After that, the central density does not 
increase securely and keeps its values p c ~ 10~ 8 g cm~ 3 , although it increases temporally 
by non-axisymmetric deformation of the first core. Note that the mass of the first core 
increases steadily by mass accretion. It seems that such non-axisymmetric instabilities in 
the first core lead to the re-distribution of the angular momentum to support larger masses 
[we will discuss this separately in Saigo, Matsumoto, & Tomisaka (2005)]. We expect that 
this massive first core fragments finally when it reaches the true maximum (mass) beyond 
which no hydrostatic configuration exists. A large-scale three-dimensional calculation must 
be performed to understand the evolution in later stages. 

On the other hand, the rotating cloud core with £l c tg > 0.03 at the central density 
of p c = 1.0 x 10~ 19 g cm~ 3 fragments promptly in the early stage of the first core phase 
(MH03). Here, we refer to their High- Resolution- Model with D, c ts = 0.046 at the central 
density of p c = 1.0 x 10 _19 g cm~ 3 as a typical fragmentation model (hereafter we refer to it as 
'HRM005'). This rapidly rotating cloud forms the first core with a mass of M core ~ 0.02M© 
at p c ~ 6 x 10~ n g cm~ 3 (Fig.76 of MH03). We plot it in our Figure 7 with a filled 
circle symbol. The position is in a non-axisymmetric unstable region on the evolution path 
with uj ~ 0.06. Soon after, this first core deforms into a bar and fragments into a binary 
(Fig. 7c? of MH03). These fragments are self-gravitating and these high-density clumps are 
embedded in the extended first core. To follow and understand the subsequent evolution 
of the fragments in our Figure 7, we define the mass of the clump, M c i ump , integrating 
the gas with p > 3.8 x 10~ n g cm~ 3 . The bar fragments into two clumps with masses of 
Afdump = 1.0xl0" 2 M o and l.lxlO~ 2 M Q atpc ~ 4xl0~ 10 g cm . Each fragment accretes 
mass and grows M c i ump = 3.0 x 1O~ 2 M and 2.4 x 1O~ 2 M at p c ~ 4 x 10~ 8 g cm" 3 (Fig. 7/ 
of MH03). We plot these fragments in Figure 7 with open circles. This evolution can be 
understood as follows: in the fragmentation process the specific angular momentum for spin 
motion from each fragment is much reduced compared with the orbital angular momentum 
and as a result, the central density of each clump rises rapidly. After a while, each fragment 
enters the non-axisymmetric unstable region again by mass accretion. It is consistent with 
the appearance of spiral arms in each fragment shown in Figure. 7f of MH03. 

It seems that through the fragmentation process a significant part of the angular momen- 
tum of the first core goes into the orbit angular momentum and the spin angular momentum 
of fragments is reduced. To estimate the angular momentum redistribution in the fragmen- 
tation process, we introduce a model core with power-law distribution of the surface density, 
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E(r) and rotation velocity v^{r) as 



S(r) oc r a , 
v^(r) oc r&, 



(12) 



where r denotes the cylindrical radius. And we assume that the gas in the disk in an opening 
angle of 0f forms a fragment without any tidal interaction. The ratio of the spin-specific 
angular momentum of fragment, jf, sp in/ jcore, becomes 



(see Appendix B). Here, we consider the first core with the uniformly rotating and uniform 
surface density (a = 0, and (5 — 1) as a model core. If the whole disk fragments into equal 
mass binary (0f = n), the ratio jf iSp in/jcorc becomes 0.64. The ratio decreases to 0.125 for 
0f <S 1. Thus the spin-specific angular momentum, j, decreases to 10% - 60% of the specific 
angular momentum of the first core before fragmentation jf !Sp in- This value is consistent with 
the results of three-dimensional simulation. For example, Machida et al. (2005) obtained 
the ratio of jf, sp in /jcore = 0.1 — 0.5 from their three-dimensional simulations of the collapse 
of six magnetized clouds (the fragment is defined as high density clumps with n > 2.5 x 10 11 



Finally, we discuss the possibility that the fragmentation induces the second collapse. 
We take the fragmentation model of MH03 shown before as the initial first core (circle 
symbol in Fig. 7). From equation (6), uo of the fragments changes in proportion to the 
change in specific angular momentum, jf )Sp i n / jcore and in inverse proportion to the fraction 
of mass contained in fa . Therefore, after the fragmentation, each self-gravitating fragment 
forms a new rotating first core with rotation speed of (0f/27r) 2 (jf iS p in /j corc )c<j. We consider 
the lowest limit of oo as the value of model RHM003 which is a no-fragmentation model 
of MH03. If it fragments into an equal mass binary, if jf )Sp i n < 0.05 j cor e, the fragments 
begin the second collapse. This redistribution fraction is smaller than the usual value of 
jf.spin = (0.1 - 0.5) j core . 

As a result, although the fragmentation can reduce the spin angular momentum of the 
cores, it is hard to induce the second collapse in the early stage of the first core. 




= 1 - 



2 sin(0 f /2) 1 2 (2 + a) (3 + a + g) 
fa J (3 + a) (2 + a + f) 



(13) 



cm 3 ). 



5.3. Effect of Magnetic Breaking 



Basu & Mouschovias (1994) pointed out that magnetic braking plays an important role 
in the case of the magnetized cloud. Tomisaka (1998, 2000, 2002) followed the collapse of 
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rotating magnetized cloud cores with axisymmetric magnetohydrodynamical (MHD) simu- 
lations using a two-dimensional nested grid method. At the end of the isothermal runaway 
collapse, the central region within ~ 0.01M Q has a rotation rate of u ~ 0.071 (Fig.l of 
Tomisaka 2000). During the first core accretion phase, the azimuthal magnetic tension 
transfers its angular momentum from the disk mid-plane to the surface and gas with excess 
angular momentum near the surface being ejected by the centrifugal force. The outflow is 
accelerated by the magnetic tension and pressure. This shows that the magnetized rotating 
first core decreases its angular momentum and continues to contract till the second collapse 
begins as well as the non-rotating spherical cloud. Tomisaka (2002) mentioned that about 
99% of the angular momentum is removed by the outflow with 10% in mass during the first 
core phase and the first core begins the second collapse. Thus, the magnetic braking moves 
the equilibrium position on the p c0 — M corc plane to the right. 

It should be noted that his calculations assumed an ideal MHD, although the coupling 
between the magnetic field and the gas weakens over the density p > 4 x 10~ 12 g cm~ 3 
(Nakano & Umebayashi 1986). Because the ionization fraction of high-density gas is quite 
low, the electric conductivity decreases as collapse proceeds. In the late phase of the first 
core, the magnetic braking may become ineffective. Therefore, the ideal MHD calculation 
may overestimate the effect of the magnetic field. On the other hand, non-magnetic hydro- 
dynamical simulations give the opposite limit. 

5.4. New Evolution Stages Induced by Rotation 

5.4-1- Another Evolutionary Path to Form the 2nd Core without Second Collapse Phase 

To begin the second collapse, it is necessary that in the first core a significant amount of 
angular momentum is extracted at once from the central region. Here, we propose that there 
are other possibilities for formation of a stellar core without the second runaway collapse 
phase. There are two kinds of evolution. First, the central density reaches the stellar density 
while the first core continually loses angular momentum due to fragmentation or gravitational 
torque. Second, the gas with small angular momentum accretes on the first core in the later 
stage and the first core contracts gradually until the stellar core is formed. 

5-4-2. The Centrifugally Supported Disk 

In our calculation, the dynamically collapsing first core forms a rotationally supported 
1.5th core before formation of the second core in models with 0.012 < u < 0.015 at p c = 
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4p diss . Bonnell & Bate (1994b) found the centrifugal bounce of the second core at the end of 
the second collapse phase in the models far from virial equilibrium initially (small rotational 
support). This leads to a self-gravitating ring with p < p stc i- The centrifugally supported 
1.5th core may lead to a binary formation and may change the mass accretion rate onto the 
second core (or star) in the early phase. 



5.5. Dependence of Angular Momentum Distribution 

We considered the evolution of the first core by assuming the angular momentum distri- 
bution of equation (5). Here, we evaluate the effect of the angular momentum distribution. 
Figure 8 shows the p c o — M core diagram for the cloud with angular momentum of a uniformly 
rotating, uniform-density cylinder, 

mm = (w) 



-^core / -^core 



Comparing Figure 7, there are no essential differences in evolution path, gravitational 
stability, and the distribution of Therefore, the mass-central density relation is 

insensitive to the angular momentum distribution. 



5.6. The Luminosity Evolution of Rotating First Core 

The evolution of luminosity and the time-scale of the first core are important quantities 
for the observation. Figure 9 shows the evolutionary curves of luminosity for models of 
u = 0.0,0.012,0.015,0.016,0.02,0.03, and 0.04. The luminosity of first core L corc is obtained 
by an energy change as 

L corc = -j t [W + E Kin + E T ], (15) 

where W, £kin and Et denote the gravitational energy, kinematical energy and thermal 
energy of the first core, respectively. Here, we assume the constant mass accretion rate of 
M = 1.0 x 10~ 5 M /yr as a reference. It is the typical mass accretion rate of a runaway 
collapsing rotating disk (Saigo & Hanawa 1998). The time t from the formation epoch 
is t — M COTe /M. Therefore, the t and L corc can be scaled by the mass accretion rate, 
t oc yt~ an d -^core oc M core . In the non-rotating case of uj = 0, the time-scale of the 

M core 

first core contraction phase is about l,000(M CO rc/A^rcf) _1 yr- The luminosity of the first core 
increases up to 0.03 (M core /M ref ) L at the end of the first core phase. This is consistent 
with (Masunaga, Miyama, & Inutsuka 1998). 
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In the case of a rotating first core, the time- scale becomes much longer. However, its 
luminosity becomes smaller for larger uo, because the central density becomes smaller for 
larger uo and the gravitational potential becomes shallower for larger uo. The luminosity 
of the rotating model of 00 — 0.03 is L CW(i ~ 0.01 (M corc /M rc f) L . On the other hand, 
the luminosity curves for models with uo = 0.012 and 0.015 show the rapid and temporal 
increase at t ~ 1500 — 2000 yr in Figure 9. This corresponds to the rapid contraction around 
p c0 ~ p 4 diss (see the horizontal path of fig. 7). 

The usual first core has a rotation with uo > 0.015. So actually the time-scale of the 
first core is several thousand years or more depending on the accretion rate. The time-scale 
of the prestellar core of nH2 ~ 10 5 cm -3 is ~ 0.1M yr (Onishi et al. 1998). Therefore, we 
expect that at least several percent of prestellar cores contain first cores and this percentage 
is inversely proportional to the accretion rate. During such a long first core phase, they must 
suffer from non-axisymmetric instability and this causes rapid contraction of first cores (see 
evolution path of MRN003 in fig. 7). Therefore, we expect that first cores repeat temporal 
increases in luminosity during their long first core phase, though it has a small average 
luminosity. 

6. SUMMARY 

We investigated the evolution of the rotating cloud during the first core and the second 
collapse phases. The evolution path of first core was constructed as a sequence of the 
equilibrium solutions. The evolution is characterized by the rotation rate uo. A cloud with 
an extremely slow rotation of uo < 0.015 evolves similarly to the spherical collapse during 
the first core phase. On the other hand, a cloud with uo > 0.015 has a certain maximum 
density and the first core does not begin the second collapse directly. The first core grows 
more and more in mass by mass accretion and gradually becomes non-axisymmetrically 
unstable. Our evolution path agrees well with the results of three-dimensional calculations. 
After that, the evolution of the first core is characterized by non-axisymmetric evolution, 
such as a formation of massive spiral arms, deformation into a bar, or fragmentation. Our 
p c o — M corc diagram helps to explain the evolution comprehensively. We expect that the 
first core with rotation does not begin the second collapse phase without loss of a significant 
amount of angular momentum. Therefore, there is another possibility that the first core 
contraction phase continues gradually till formation of the stellar core. During the second 
collapse phase, a cloud with 0.015 < uo < 0.012 forms a centrifugally supported disk before 
the stellar core formation. 

We expect that the rotating first core has a the time-scale of several thousand years or 
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more and at least several percent of prestellar cores contain first cores. The rotating first 
core has a small average luminosity of L corc = 0.003 — 0.03(M core /M re f) L . However, we 
expect that the rotating first core repeats temporal increases in luminosity during its first 
core phase. 

We wish to express our gratitude to I. Hachisu who offered his numerical code (SCF) 
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A. RELATION BETWEEN NON-DIMENSTIONAL ANGULAR SPEED 

AND THE ROTATION RATE oo 

The rotation rate u can be rewritten using the central angular speed fl c and central 
density p c , when the disk has the uniform rotation, uniform surface density, E and isothermal 
hydrostatic equilibrium in the ^-direction. In this case, the rotation rate is 

\^2,C[ so t/ C ore V^Qso^c / » -, \ 

where Cj SO is the isothermal sound speed. In the isothermal hydrostatic equilibrium, the 
uniform surface density, E, is related to the central density as 

E = (^l^:)V2 (A2) 

7rG 

If we define the free-fall time-scale as t s = l/(47rGp c ) 1 / 2 , the rotation rate is rewritten by 
non-dimensional angular speed, 

uj = n c t s . (A3) 
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B. DISTRIBUTION OF ANGULAR MOMEMTUM TO ORBITAL AND 

SPIN 

We assume a simple model of the first core with power-law distributions of the surface 
density, S(r) and rotation velocity, v^r) as 

S(r)=Ar Q , 

Vip (r) = Brf 3 , { ' 

where r denotes the cylindrical radius. For a > — 2, mass of the disk is expressed as 

/■■Rcoro ITT A 

M corc = / 2vrrS(r)rfr = ——H£*, (B2) 
Jo ol + 2 

where -R cor e is the outer radius of the first core. The average specific angular momentum of 
the first core becomes 

Jcorc " m^ c " ^TJT^ BR -- (B3) 

We assume that the disk fragments and gas in the opening angle of 0f forms a compact object 
(fragment) without any tidal interaction. The mass, position and velocity of the fragment 
are M f = (0 f /27r)M core , r f = J corc S(r) r rdrd<f)/M u and v t = J corc S(r-) v rdrd<f)/M u 
respectively, where we take the foot of the position vector to be at the disk center. They 
can be rewritten as 



rRcore r(n+^)/2 2 sin(0f/2) a + 2 „ 

r t = r { = Mf 1 / / E r r r sin <f>dj>dr = V ; ■ A rc 

Jo J(?-to)/2 <Pf « + 3 

i f Rc ° IC /" (7r+ ^ )/2 w n / «. . JfJJ 2sin(0 f /2) a + 2 , 

Jo Ju-6 f )/2 ' <Pt a + [3 + 2 



v\\ = 

(B4) 

where v± and v\\ denote the perpendicular and parallel components of velocity to the position 
vector r, respectively. The specific orbital angular momentum of the fragment is jf )0r bit — 
r x v±. Therefore, the ratio of the spin-specific angular momentum of fragment, jf, sp in/jcore, 
becomes 



Ji. 



spm 



1 - 



2sin(0 f /2) 

0f 



(2 + a)(3 + a + /3) 

(3 + a)(2 + a + (3)- 1 ' 
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r [AJ| } [AUf > |AUJ 

Fig. 1. — Rotating first cores with the same central density of p C Q = 4pdiss but different rota- 
tion speeds. The contours and gray scale denote the density distribution in logarithmic scale 
on the r-z plane. The white dashed curve denotes the isodensity contour of constant density 
with dissociation density pdiss- They are examples for (a) slow rotation with T/|W| = 0.11 
(Model A), (b) moderate rotation with = 0.15 (Model B), and (c) fast rotation with 

T/lV^I = 0.22 (Model C), where T"/|W| is the ratio of rotational energy to the gravitational 
energy. 
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10" 11 10" 10 10" 9 10" 8 10" 7 10" 6 10~ 5 
PcO ( 9 cm " 3 ) 

Fig. 2. — The central density core mass (p($-M COTe ) diagram. The evolution path is plotted 
as sequences of equilibrium solutions. The solid curves denote the evolution path and each 
path is labeled by rotation rate u. The thick curve denotes the evolution path of the non- 
rotation spherical model. The red region shows the gravitationally unstable region. Circles 
denote first cores which form second cores directly. Triangles denote first cores which form 
centrifugally supported disks before the formation of the stellar core. Crosses indicate stable 
models with no contraction. 
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Fig. 3. — Density and velocity distributions on the r-z plane for Model A at six different 
stages. The central density is (a) p c = 2.3 x 10~ 7 g cm" 3 , (6) 2.2 x 10" 6 g cm -3 , (c) 2.9 x 10~ 5 
g cm -3 , (d) 2.8 x 10~ 4 g cm -3 , (e) 1.81 x 10~ 3 g cm" 3 , and (J) 1.0 x 10~ 3 g cm" 3 at each 
stage. The contour curves and grayness denote the density in the logarithmic scale. The 
arrows show the velocity distributions. 
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Fig. 4. — (a)density, (b) radial infall velocity, and (c)rotation velocity on the equatorial plane 
are shown for model A. Panel (d) shows the infall velocity along the z- axis. Solid curves 
represent the distribution in the second collapse phase. The evolution after the bounce is 
shown with dashed curves. Curves (a)-(e) correspond to the snapshots of Fig. 3 




Fig. 5. — Same as Fig. 3 but for model B. The white dotted curves denote the 1.5th core. 
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Fig. 6. — Same as Fig. 4 but for model B. Curves (a)-(a) correspond to the snapshots of 
Fig. 5. 
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10" 11 10" 10 10" 9 10" 8 10" 7 10" 6 10" 5 
PcO ( 9 cm " 3 ) 

Fig. 7. — Same as Fig. 2 but with rotational energy and three-dimensional hydro results. 
Color denotes the ratio of the rotational energy to the gravitational energy, T/|W|, and 
the red dashed curve denotes the critical value of T/|VV| = 0.273 for the non-axisymmetric 
instability. The squares denote the results of Bate (1998). The filled and open circles denote 
the results of HRM005 of MH03. The first core for HRM005 model fragments into binary 
(open circles). The white dashed curve denotes the results for HRM003 of MH03. 
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Fig. 8. — Same as Fig. 7 but for a different angular momentum distribution. The angular 
momentum distribution for the uniformly rotating, uniform density cylinder is assumed. 
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Fig. 9. — Evolutionary curve of the first core luminosity during the first core phase. A 
constant mass accretion rate of M — 1.0 x 10~ 5 M /yr is assumed as a reference. 



